Patterns of interval correlations in neural oscillators with adaptation 
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Neural firing is often subject to negative feedback by adaptation currents. These currents can 
induce strong correlations among the time intervals between spikes. Here we study analytically the 
interval correlations of a broad class of noisy neural oscillators with spike-triggered adaptation. Our 
weak-noise theory provides a general relation between the correlations and the phase-response curve 
of the oscillator, gives intuition about the sign of the correlations, and predicts a variety of correlation 
patterns that depend on the nonlinear dynamics of the model. At high firing rates, the long-term 
variability associated with the cumulative interval correlations becomes small, independent of model 
details. Our results are verified by comparison with stochastic simulations of the exponential, leaky, 
and generalized integrate-and-fire models with adaptation. 



The nerve cells of the brain generate action potentials 
(spikes) by a nonlinear, adaptive, and noisy mechanism. 
In order to understand neural signal processing in single 
neurons it is vital to analyze the sequence of the inter- 
spike intervals (ISIs) between adjacent action potentials. 
There is experimental evidence accumulating that the 
spiking in many cases is not a renewal process, i.e. a 
spike train with mutually independent ISIs but that in- 
tervals are typically correlated over a few lags [HE]- Such 
correlations are often characterized by the serial correla- 
tion coefficient (SCC) 

^ ((T,-(T,))(r,+,-(r,+,))) 

where Ti and Ti+k are two ISIs lagged by an integer k 
and (•) denotes ensemble averaging. 

An ubiquitous mechanism for ISI correlations are slow 
feedback processes mediating spike-frequency adaptation 
[31 13] . These adaptation currents are typically associated 
with short-range correlations with a negative SCC at lag 
k — 1 and a reduced Fano factor as demonstrated by sev- 
eral numerical [SHS] and analytical studies [9lfl2]. How- 
ever, a simple theory that predicts and explains possible 
correlation patterns is still lacking. 

In this letter, we present a relation between the ISI cor- 
relation coefficient pk and a basic characteristics of non- 
linear neural dynamics, the phase-response curve (PRC). 
The PRC quantifies the advance (or delay) of the next 
spike caused by a small depolarizing current applied at 
the time t after the last spike. For neurons which in- 
tegrate up their input (integrator neurons), the PRC is 
positive at all times (type I resetting) whereas neurons, 
which show subthreshold resonances (resonator neuron), 
possess a PRC that is partly negative (type II resetting) 
[13] . Below we show that resonator neurons possess 
a richer repertoire of correlation patterns than integrator 
neurons do. 

Model. Spike frequency adaptation can be modeled 
by Hodgkin-Huxley type neurons with a depolarization- 
activated adaptation current [H [TSl [Hj. However, the 



spiking of such conductance-based models can in many 
instances be approximated by simpler multi-dimensional 
integrate-fire (IF) models that are equipped with a spike- 
triggered adaptation current [E]; adapting IF models 
perform excellently in predicting spike times of real 
cells under noisy stimulation |18| . Here, we consider 
a stochastic nonlinear multi-dimensional IF model for 
the membrane potential v, N auxiliary variables Wj 
{j = 1,...,N) and a spike-triggered adaption current 
a{t): 

V = fa{v,w) + 11- a + ^(t), Wj = fj{v,w), (2a) 
Taa= -a + TaAY,S{t-U)- (2b) 

i 

The membrane potential v{t) is subject to weak Gaus- 
sian noise ^{t) with (C(t)^(i')) = 2D6{t - t') and noise 
intensity D. The dynamics is complemented by a spike- 
and- reset mechanism: whenever v{t) reaches a threshold 
wt, a spike is registered at time U = t and v(t) and 
w(<) = [wi{t), . . . ,WN{t)]'^ are reset to v{ti+) = and 
w(<i+) — w^; a{t) suffers a jump by A as seen from 
Eq. ( |2b| , which resembles high-threshold adaptation cur- 
rents^! [3 • The constant input current /z is assumed to 
be sufficiently large to ensure ongoing spiking even in the 
absence of noise. 

An important special case, the adaptive exponential 
integrate-and-fire model is illustrated in Fig.l (here, 
/o(w) = -^v + 7ATexp[(u - 1)/At] and N = 0). Time 
courses of v{t) and a{t) are shown in Fig.lal,bl for two 
distinct correlation patterns possible in this model. The 
ISI sequence Ti,Ti+i,Ti+2 displays patterns of short- 
long-long (Fig.lal) and short-long- short (Fig.lbl), cor- 
responding to a negative SCC, which decays monotoni- 
cally with the lag k (Fig.laS) or to an SCC oscillating 
with k (Fig.lbS). Below we develop a theory to analyze 
these and other correlation patterns possible in multi- 
dimensional adapting IF models. 

Theory. In our model Eq. ([2]), a{t) is the only variable 
that keeps a memory of the previous spike times thereby 
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FIG. 1: Correlation patterns in the adaptive exponential IF 
model with t„ = 10,7 = 1,At = 0.1, ut = 2, D = 0.1. 
Adaptation is weak (A = 1,/i = 15) in (a) and strong 
(A = 10, = 80) in (b). Membrane voltage v{t) and adap- 
tation variable a(t) with ISI sequences {Ti} and peak adap- 
tation values {ai} are shown in (al,bl). Colored pieces of 
trajectories in the phase plane {v, a) in (a2,b2) correspond 
to the respective colors in (al,bl); deterministic limit cycle 
(LC) is indicated by a thick black line. For weak adaptation 
(a2) a short ISI Ti causes positive deviations 5ai — at — a* 
and Sui+i = a^+i — a* of peak values leading to long ISIs 
Ti+i and Tn.2 and, hence, to a negative ISI correlation at all 
lags (a3). Because of the qualitatively different limit cycle 
for strong adaptation (b2), deviations Sui and (Ja^+i differ in 
sign, yielding an oscillatory correlation pattern (b3). 



inducing correlations between ISIs. Over one ISI the time 
course of adaptation is an exponential decay, relating two 
adjacent peak values ai and aj+iof a(t) by 



a,:e 



A. 



(3) 



We assume that in the deterministic case {D — 0) our 
model has a finite period T* (i.e. the model operates in 
the tonically firing regime) and, hence, for D = the 
map ([3| has a stable fixed point 



A/ [1 - exp(-T7Tj] 



(4) 



The asymptotic deterministic dynamics can be inter- 
preted as a limit-cycle like motion in the phase space 
from the reset point to the threshold and back by the 
instantaneous reset (cf. Fig.la2,b2). 

Weak noise will cause small deviations in the period 
6T, =T,-T*mT,- (Ti) that are mutually correlated 
with coefficient pk- The peak adaptation values, however, 
also fluctuate, 6ai — ai — a*, and both deviations are 
related by linearizing Eq. ([3|: 



a* \ 



(5) 



A second relation between the small deviations can be 
gained by considering how a small perturbation affects 
the length of the period. This effect is captured by the in- 
finitesimal phase response curve (PRC), Z{t), t G (0,T*) 
[131 [13]. In our case, the perturbation ^(i) — Jai exp(— (t— 
ti)/Ta) comprises the weak noise and the deviation in the 
adaptation during (i,, ij+i): 

6T,+i^J^ dtZ{t)(^5a,e-^ -^{U + t)y (6) 

Combining Eqs. ([s]), ([6]) we obtain the stochastic map 

Sa,+i = a-dda^ + E,, (7) 

where = —^^^ dtZ(t)^{ti + t) are uncorrelated 
random numbers and 



'a Jo 



dtZ{t)e~ 



(8) 



Note that local stability of the fixed point a* requires 
that |q!??| < 1. The covariance Ck = (daiSai^k) of the 
auto-regressive process Eq. ([7| can be calculated by ele- 
mentary means and using Eq. ^ we obtain for k > 1: 



p^ = -A{l-^){ad) 



k-l 



A 



1 + a2 - 2aH' 



(9) 



In order to compute a and d via Eq. ([s]), we have to 
calculate T* and Z{t) (a* then follows from Eq. Q), 
which can be done for simple systems analytically. 

Our main result, Eqs. ([8]),(|9]), shows that the SCC is 
always a geometric sequence with respect to the lag k 
that can generate qualitatively different correlation pat- 
terns depending on the value of ■& and thus on PRC and 
adaptation current. Because \ad\ < 1 and < a < 1, the 
prefactor A in Eq. ^ is always positive. Consequently, 
pi is negative for i? < 1 and positive for {) > 1. The sign 
of higher lags is determined by the base of the power: 
for z9 > correlations decay monotonically, whereas for 
I? < the SCC oscillates. Two special cases are d — 
with a negative correlation at lag 1 and vanishing corre- 
lations at all higher lags and ■& — 1 where all correlations 
vanish. Our geometric formula covers all patterns of in- 
terval correlations that have been theoretically studied 
[H [19] , including the perfect IF model plj . 

The cumulative effect of the correlations can be de- 
scribed by the sum over all pk that is related to the long- 
time limit of the Fano factor and the low-frequency limit 
of the spike train power spectrum. It is given by 
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(10) 



where the latter expression is valid at high firing rates 
with T* <ti Ta, achieved by a strong input current p. In 
particular, for strong adaptation (Atq ^ vt) the sum is 
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only slightly larger than —1/2 (for —1/2 the Fano factor, 
and hence the long-term variability vanishes). 

One- dimensional IF models with adaptation. In the 
simplest case {N — 0, /o(w,w) = f{v)) the PRC reads 

Z{t)^Z{T*)exp[£ dt' fivoit'))], (11) 

where vo{t) is the limit cycle solution and Z{T*) = 
[/(ut) + h — a* + A]~^ is the inverse of the velocity 
Vo{T*) at the threshold, which is always positive. Thus, 
the PRC is positive for all t e (0,r*), i.e. ID-IF models 
show type I behavior. Looking at Eq. Q we find that a 
positive PRC inevitably yields -i^ < 1 and hence implies 
pi < 0. Intuitively, a short ISI causes in the following on 
average a higher inhibitory adaptation during the subse- 
quent ISI. Such an inhibitory current always enlarges the 
ISI in type I neurons - hence, a short ISI is followed by 
a long ISI. 

The sign of the correlations at higher lags can be in- 
ferred from the sign of i?, for which one can show that 
■ff = (/(O) + ^ - a*)Z{0). Because Z{0) > 0, the sign of 
i9 is determined by the sign of /(O) + fi — a* . For suf- 
ficiently small A (weak adaptation, see Fig. [l^) we will 
have a* < /(O) -I- consequently, > and a negative 
correlation at all lags (Fig. [1^3). In this case, a short ISI 
occurring by fluctuation will cause a positive deviation 
6ai (Fig. [1^2, green arrow). Geometrically, it is plausible 
that such a positive deviation causes a likewise positive 
deviation 5ai+i in the subsequent cycle (Fig. [1^2, red ar- 
row). Because a positive deviation is associated with a 
long ISI, the initial short ISI is on average followed by 
longer ISIs. 

In marked contrast, if A is sufficiently large (strong 
adaptation) such that a* > /(O) + fJ,, "d becomes negative 
and hence the SCC's sign alternates with the lag. This 
alternation of the sign can be understood by means of 
the phase plane. Let us again consider a positive devi- 
ation Soi due to a short preceding ISI (Fig. [i|d2, green 
arrow). Because i'o(O) — /(O) + fi — a* < 0, the neuron 
hyperpolarizes at the beginning of the interval (i.e. the 
trajectory moves left into the region of negative voltage). 
At the turning point of the limit cycle the deviation of 
the adaptation reverses its sign, and hence Ja^+i becomes 
negative (Fig.[l|^b2), red arrow). Because a positive (neg- 
ative) deviation corresponds on average to a long (short) 
ISI, the alternation of Soi also entails an alternation of 
the ISI correlations. 

As demonstrated in Fig{T^3,b3, our theory works well 
for the adapting exponential integrate-and-fire model. 
We next demonstrate the validity of our approach over a 
broad range of firing rates (Fig. [2]) for another important 
ID model, the adapting leaky integrate-and-fire model 
[20] for which f{v) = —jv and 

Z{t) = exp[7(< - T*)]/{fi - jVT ~a*+A) (12) 
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FIG. 2: ISI correlations of the adapting LIF model vs. firing 
rate l/(Ti) ~ 1/T* where the rate is varied by increasing ^. 
The panels display (from top to bottom) pi, p2 and the sum 
X^fcLi Pk for simulation (circles, m — 100) and theory (solid 
lines, m — >■ oo). (a) Moderate adaptation: A = 1, (b) strong 
adaptation: A = 10. Both: 7 = 1, Ta = 10, D = 0.1, vt = 1. 



(here T* has still to be determined from a transcenden- 
tal equation). Changing the firing rate by varying the 
input current p, we find a good agreement for the first 
two correlation coefficients and the sum of all pk- In 
accordance with previous findings [TTJ [111 HI] j the 
first correlation coefficient pi displays a minimum cor- 
responding to strong anti-correlations between adjacent 
intervals. The correlations at lag 2 can be positive for 
a finite range of firing rates if the adaptation strength 
A is sufficiently large (Fig. [2](b)), whereas for moderate 
A we find a negative p2 at all firing rates (Fig. [2][a)). 
In both cases, however, the sum of SCCs approaches a 
value close to —1/2 for high firing rates as predicted by 
Eq. ( 10 1 (Fig.[2j bottom). This is strikingly similar to ex- 
perimental data from weakly electric fish, in which some 
electro-receptors display a monotonically decaying SCC 
and some show an oscillatory SCC j4« but all cells exhibit 
a sum close to —1/2 [22) . 

Adapting GIF model. Different correlation patterns be- 
come possible if we consider a type II PRC, which is 
by definition partly negative and can lead to a negative 
value of the integral in Eq. (|8]), and hence to d > 1. 
This corresponds to a non-negative SCC at lag 1, which 
is infeasible in the one-dimensional case. To test the pre- 
diction pi > 0, we study the generalized integrate-and- 
fire (GIF) model [33] (also known as resonate-and-fire 
model |;24 ), which is defined by fo{v,w) = ~jv — /3w 
and fi{v,w) = v — w. For this model, we find 
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where 1/ = 7 l/r^, O = 
one component of the deterministic limit-cycle solution 



(13) 

^ and wo{t) is 
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[vo{t),wo{t),ao{t)] that we calculated numerically. 

(a) integrator (LIF) (b) resonator (GIF) 
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FIG. 3: Possible correlation patterns and corresponding 
PRCs. For the adapting LIF model (a), i? < 1 and only three 
qualitative difterent cases are possible. The adapting GIF 
model (b) exhibits the full repertoire of correlation patterns 
because the PRC can be partly negative and i? can attain val- 
ues from its entire physically meaningful interval [— 1/a, l/a]. 
The value of i? and hence the type of correlation pattern is set 
by the integral over the weighted PRC Z{t) = Z(t)e"^ 
shown in left panels. LIF parameters: D = 0.1, Ta = 2, (i) 
M = 20, A = 10, (ii) M = 20, A = 4.47, (ih) ^ = 5, A = 1. 
GIF parameters: (i) = 10, /3 = 3, = 10. (u) = 11.75, 
13 = 3, ra = 10. (iii) M = 20, /3 = 1.5, r, = 10. (iv) = 2.12, 
13 = 1.5, Ta = 1, A = 10. (v) At = 1.5, /3 = 1.5, t„ = 1, A = 9 
D = 10~^. Unless stated otherwise, 7 = 1, A = 1, r^, = 1.5, 

D = 10"", Wr = 0. 

In Fig. [3]d we demonstrate that all possible correla- 
tion patterns can be realized in the GIF model and that 
the predicted SCCs agree quantitatively well in theory 
and model simulations (for comparison, see the SCC for 
the LIF in Fig. [3^). To each distinct pattern belongs 
a range of (Fig. |3j left), determined by the area un- 
der the weighted PRC Z{t) = '^e~^Z{t). The func- 
tion Z{t) (left column in Fig. |3^,b) illustrates, why an 
adapting GIF neuron can show vanishing (Fig. [3]3(iv)) 
or even purely positive ISI correlations (Fig. [3}d(v)). In 
case of type II resetting, inhibitory input can shorten the 
ISI because of the negative part in the PRC; here in- 
hibition acts like an excitatory input. Consequently, a 
short ISI will induce a stronger inhibition (adaptation) 
that now causes a likewise short interval and results thus 
in a positive correlation between adjacent ISIs. Also, 
the shortening effect of the adaption current in the early 
negative phase of the PRC can be exactly balanced by 
the delaying effect of the late positive phase of the PRC 



(pseudo-renewal case, in which the area under Z is zero). 

Conclusions. We have found a general relation between 
two experimentally accessible characteristics: the serial 
interval correlations and the phase response curve of a 
noisy neuron with spike-triggered adaptation. Our ana- 
lytical results require that the noise is weak but are valid 
for arbitrary adaptation strength and time scale. The 
theory predicts distinct correlation patterns like short- 
range negative and oscillatory correlations that have been 
observed in experiments [H and in simulation stud- 
ies of adapting neurons [H [7] . Beyond negative and os- 
cillatory correlations, we have found, however, that res- 
onator neurons with spike-frequency adaptation can ex- 
hibit purely positive ISI correlations or a pseudo-renewal 
process with uncorrelated intervals. Adaptation currents 
that are commonly associated with negative ISI correla- 
tions can thus induce a richer repertoire of correlation 
patterns than previously thought. Despite the multitude 
of patterns, there is a universal limit for the cumula- 



tive correlations at high firing rates (cf. Eq. (lOl), which 
shows that the long-term variability is in this limit always 
reduced in agreement with experimental studies [22j . 

As an outlook we sketch, how our theory could be used 
to constrain unknown physiological parameters by mea- 
sured correlation coefficients and phase response curves. 
For instance, from the mean ISI we can estimate T* — 
(T). Furthermore, knowing pi = —A(a, i?)(l — "i?) as well 
as the ratio P2/P1 = (^"^ one can eliminate i? and solve 
for a. This allows to estimate the unknown adaptation 
time constant Tq — —T*/\iia and the amplitude of the 
adaptation current 




dtZ{t)e~^ 



(14) 



Although experimental PRCs are notoriously noisy |13j, 
the integral over Z(t) determining our estimate of a* is 
less error-prone. Combining our approach with advanced 
estimation methods for the PRC f26^, may thus provide 
an alternative access to hidden physiological parameters 
using only spike time statistics. 
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